function pairedGroup = autoJoinPZ(polos,zeros)

    polosDob = polos((imag(polos)>0));
    polosReal=polos(imag(polos)==0);
    zerosDob = zeros((imag(zeros)~=0));
    zerosReal = length(zeros(abs(zeros) == 0));
    pairedGroup = [];
    
    pairGroup = struct('polo1','polo2','zero1','zero2');
    while ~isempty(polosDob)
        dist = inf;
        index = getMaxQIndex(polosDob);
        zeroindex = 0;
        withReal = 0;
        
        for k = 1:length(zerosDob)
            if( abs(polosDob(index)-zerosDob(k)) < dist)
                dist = abs(polosDob(index)-zerosDob(k));
                zeroindex = k;
            end
        end
        
        if abs(polosDob) < dist

            if zeroindex == 0 && zerosReal == 1
                zerosReal = 0;
                withReal = 1;
            end
            
            if zerosReal >= 2
                zerosReal = zerosReal - 2;
                withReal = 2;
            end
        end
        
        pairGroup.polo1 = polosDob(index);
        pairGroup.polo2 = conj(pairGroup.polo1);
        polosDob(index) = [];
        
        if withReal == 1
            pairGroup.zero1 = 0;
            pairGroup.zero2 = [];
        else
            if withReal == 2;
                pairGroup.zero1 = 0;
                pairGroup.zero2 = 0;
            else
                if zeroindex ~= 0
                    pairGroup.zero1 = zerosDob(zeroindex);
                    pairGroup.zero2 = conj(pairGroup.zero1);
                    zerosDob(zeroindex) = [];
                else
                    pairGroup.zero1 = [];
                    pairGroup.zero2 = [];
                end
            end

        end
        pairedGroup = [pairedGroup; pairGroup];
    end
    
    while ~isempty(polosReal)
        dist = inf;
        index = getMinRealIndex(polosReal);
        zeroindex = 0;
        withReal = 0;
        
        pairGroup.polo1 = polosReal(index);
        polosReal(index) = [];
        
        if ~isempty(polosReal)
            
            index = getMinRealIndex(polosReal);
            pairGroup.polo2 = polosReal(index);
            polosReal(index) = [];
            relW0 = sqrt(pairGroup.polo2*pairGroup.polo1);
            
            for k = 1:length(zerosDob)
                if( abs((relW0)-zerosDob(k)) < dist)
                    dist = abs(relW0 - zerosDob(k));
                    zeroindex = k;
                end
            end

            if abs(relW0) < dist
                
                if zeroindex == 0 && zerosReal == 1
                    zerosReal = 0;
                    withReal = 1;
                end
                
                if zerosReal >= 2
                    zerosReal = zerosReal - 2;
                    withReal = 2;
                end
            end



            if withReal == 1
                pairGroup.zero1 = 0;
                pairGroup.zero2 = [];
            else
                if withReal == 2;
                    pairGroup.zero1 = 0;
                    pairGroup.zero2 = 0;
                else
                    if zeroindex ~= 0
                        pairGroup.zero1 = zerosDob(zeroindex);
                        pairGroup.zero2 = conj(pairGroup.zero1);
                        zerosDob(zeroindex) = [];
                    else
                        pairGroup.zero1 = [];
                        pairGroup.zero2 = [];
                    end
                end

            end
        else
            if zerosReal == 1
                pairGroup.polo2 = [];
                pairGroup.zero1 = 0;
                pairGroup.zero2 = [];
            else
                pairGroup.polo2 = [];
                pairGroup.zero1 = [];
                pairGroup.zero2 = [];
            end
        end
        pairedGroup = [pairedGroup; pairGroup];

    end
    
end

function index = getMinRealIndex(polos)

    minpole = inf;
    index = 0;
    for k=1:length(polos)
        if abs(polos(k)) < minpole
            minpole = abs(polos(k));
            index = k;
        end
    end

end

function index = getMaxQIndex(polos)
    Q=0;
    index=0;
    for k=1:length(polos)
    
        if abs(polos(k))/abs(real(polos(k))) > Q
            Q=abs(polos(k))/abs(real(polos(k)));
            index=k;
        end
        
    end

end